/**
 * Moon phase calculation class
 * Adapted for PHP from Moontool for Windows (http://www.fourmilab.ch/moontoolw/) 
 * by Samir Shah (http://rayofsolaris.net)
 * Last modified August 2010
 
	This program is free software: you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation, either version 3 of the License, or
	(at your option) any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program.  If not, see <http://www.gnu.org/licenses/>.

 *
	USAGE: simply create an instance of the moonphase class, supplying
	a UNIX timestamp for when you want to determine the moon phase. The following
	variables will be created in the resulting object, which you can then access
	from your scripts:
	- phase: the terminator phase angle as a fraction of a full circle (i.e., 0 to 1)
	- illum: the illuminated fraction of the Moon, in degrees
	- age: the Moon's age in days
	- dist: the distance of the Moon from the centre of the Earth
	- angdia: the angular diameter subtended by the Moon as seen by an observer at the centre of the Earth
	- sundist: the distance to the sun in kilometres
	- sunangdia: the angular diameter subtended by the Moon as seen by an observer at the centre of the Earth
 **/
function moonphase(pdate){
    pdate=pdate/1000;
    this.epoch = 2444238.5;
    this.elonge = 278.833540;		
    this.elongp = 282.596403;		
    this.eccent = 0.016718;			
    this.sunsmax = 1.495985e8;		
    this.sunangsiz = 0.533128;
    this.mmlong = 64.975464;
    this.mmlongp = 349.383063;
    this.mlnode = 151.950429;		
    this.minc = 5.145396;			
    this.mecc = 0.054900;			
    this.mangsiz = 0.5181;			
    this.msmax = 384401;			
    this.mparallax = 0.9507;		
    this.synmonth = 29.53058868;	
    this.lunatbase = 2423436.0;		


    var self=this;
    this.fixangle=function (a) {
        return ( a - 360 * Math.floor(a / 360) );
    };
    this.rad2deg=function(r){
        return (r/Math.PI)*180;
    };
    this.deg2rad=function(d){
        return (Math.PI*d)/180;
    };
    this.kepler=function(m, ecc) {

            epsilon = Math.pow(1, -6);
            e = m = self.deg2rad(m);
            do {
                    delta = e - ecc * Math.sin(e) - m;
                    e -= delta / ( 1 - ecc * Math.cos(e) );
            } 
            while ( Math.abs(delta) > epsilon );
            return e;
    };
    this.pdate =  pdate / 86400 + 2440587.5;
    this.Day = this.pdate - this.epoch;								
    this.N = this.fixangle((360 / 365.2422) * this.Day);		
    this.M = this.fixangle(this.N + this.elonge - this.elongp);		
    this.Ec = this.kepler(this.M, this.eccent);					
    this.Ec = Math.sqrt((1 + this.eccent) / (1 - this.eccent)) * Math.tan(this.Ec / 2);
    this.Ec = 2 * this.rad2deg(Math.atan(this.Ec));						
    this.Lambdasun = this.fixangle(this.Ec + this.elongp);		
    this.F = ((1 + this.eccent * Math.cos(this.deg2rad(this.Ec))) / (1 - this.eccent * this.eccent));	
    this.SunDist = this.sunsmax / this.F;							
    this.SunAng = this.F * this.sunangsiz;							
    this.ml = this.fixangle(13.1763966 * this.Day + this.mmlong);				
    this.MM = this.fixangle(this.ml - 0.1114041 * this.Day - this.mmlongp);		
    this.MN = this.fixangle(this.mlnode - 0.0529539 * this.Day);				
    this.Ev = 1.2739 * Math.sin(this.deg2rad(2 * (this.ml - this.Lambdasun) - this.MM));		
    this.Ae = 0.1858 * Math.sin(this.deg2rad(this.M));								
    this.A3 = 0.37 * Math.sin(this.deg2rad(this.M));									
    this.MmP = this.MM + this.Ev - this.Ae - this.A3;									
    this.mEc = 6.2886 * Math.sin(this.deg2rad(this.MmP));								
    this.A4 = 0.214 * Math.sin(this.deg2rad(2 * this.MmP));							
    this.lP = this.ml + this.Ev + this.mEc - this.Ae + this.A4;								
    this.V = 0.6583 * Math.sin(this.deg2rad(2 * (this.lP - this.Lambdasun)));				
    this.lPP = this.lP + this.V;												
    this.NP = this.MN - 0.16 * Math.sin(this.deg2rad(this.M));							
    this.y = Math.sin(this.deg2rad(this.lPP - this.NP)) * Math.cos(this.deg2rad(this.minc));			
    this.x = Math.cos(this.deg2rad(this.lPP - this.NP));									
    this.Lambdamoon = this.rad2deg(Math.atan2(this.y, this.x)) + this.NP;						
    this.BetaM = this.rad2deg(Math.asin(Math.sin(this.deg2rad(this.lPP - this.NP)) * Math.sin(this.deg2rad(this.minc))));		
    this.MoonAge = this.lPP - this.Lambdasun;								
    this.MoonPhase = (1 - Math.cos(this.deg2rad(this.MoonAge))) / 2;					
    this.MoonDist = (this.msmax * (1 - this.mecc * this.mecc)) / (1 + this.mecc * Math.cos(this.deg2rad(this.MmP + this.mEc)));
    this.MoonDFrac = this.MoonDist / this.msmax;
    this.MoonAng = this.mangsiz / this.MoonDFrac;	
    this.phase = this.fixangle(this.MoonAge) / 360;					
    this.illum = this.MoonPhase;										
    this.age = this.synmonth * this.phase;							
    this.dist = this.MoonDist;										
    this.angdia = this.MoonAng;										
    this.sundist = this.SunDist;										
    this.sunangdia = this.SunAng;
}


